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Fully analytic formulas, which do not involve any numerical integration, are derived for the dis- 
cretized influence functionals of a very extensive assortment of spectral distributions. For Feynman 
integrals derived using the Trotter splitting and Strang splitting, we present general formulas for the 
discretized influence functionals in terms of proper integrals of the bath response function. When 
an analytic expression exists for the bath response function, these integrals can almost always be 
evaluated analytically. In cases where these proper integrals cannot be integrated analytically, 
numerically computing them is much faster and less error-prone than calculating the discretized 
influence functionals in the traditional way, which involves numerically calculating integrals whose 
bounds are both infinite. As an example, we present the analytic discretized influence functional for 
a bath response function of the form a(t) = ^ pje jt , which is a natural form for many spectral 
distribution functions (including the very popular Lorentz-Drude/Debye function), and for other 
spectral distribution functions it is a form that is easily obtainable by a least-squares fit. Evaluating 
our analytic formulas for this example case is much faster and more rigorous than numerically calcu- 
lating the discretized influence functional in the traditional way. In the appendix we provide analytic 
expressions for pj and Qj for a variety of spectral distribution forms, and as a second example we 
provide the analytic bath response function and analytic influence functionals for spectral distribu- 
tions of the form J(uj) oc uj e e~( tJ ' u ' a ' . The value of the analytic expression for this bath response 
function extends beyond its use for calculating Feynman integrals. We also provide open source 
MATLAB and Mathematica programs to make the results of this paper very easy to implement. 
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INTRODUCTION 

Influence functionals are used to describe the environment's influence on the reduced density operator dynamics (RpD) of open 
quantum systems (OQSs) in the Feynman integral representation[3l \4. l33l| - They have been used since 1959 J33J in this framework, 
which later became the basis for some of the first successful methods for calculating 'numerically exact' RpD in OQSs with large 
environments. The most popular model for OQSs is currently a Feynman- Vernon model (see equation [3]) with a continuous spectral 
distribution J(w) (see equation U) , in which the entire influence of the OQS's environment on its RpD is captured by the influence 
functional. 

The remainder of this paper will deal with this OQS model, for which discretized Feynman i nteg rals that use discretized influence 
functionals (DIFs) have been used extensively for calculating the RpD of various OQSs[9j, [l3|, [30j|. For this model, these DIFs are 
traditionally expressed in terms of improper integrals of J{uj) whose integration bounds are infinite in both directions|14l. Il5|. The 
accuracy of these integrals can be extremely important when using the Feynman integral to numerically calculate RpD in OQSs, and 
a converged calculation of them can be very computationally demanding in some cases. Aside from the inconvenience of having to 
check for convergence when numerically calculating these integrals, not having an analytic expression for the DIF makes it difficult to 
examine properties of it, such as its dependence on the parameters of the physical system and the dependence of its size on the size of 
the time steps in the discretization of the Feynman integral. 

In this paper we will give a general expression for the DIF, in terms of proper integrals of the bath response function a(t), instead of 
improper integrals of J{lo) whose integration bounds are infinite in both directions. As long as a(t) is represented analytically, these 
integrals are easy to evaluate analytically, or if they cannot be, they are still much easier to evaluate numerically than the integrals in 
the traditional expression for the DIF, since their integration bounds are finite and very tiny. 

One particular form for a(t) for which these proper integrals of a(t) can very easily be evaluated, leading to a fully analytic DII|J, is: 



K 



a(*)=J>e^. (^ 



This is a very important form for a(t) for at least two reasons. 

(1) Every physically relevant spectral distribution can be fitted to this form, and in fact many very versatile spectral distribution 
forms (including the very popular Lorentz-Drude/Debye function) naturally have bath response functions of this form (see table [J) . 

(2) The main purpose of Feynman integrals for RpD in OQSs is to benchmark methods that are less computationally expensive. 
When benchmarking, it is ideal for both methods to use the exact same form for a(t), and many of the techniques that one may wish 
to benchmark actually require a(t ) to be in this form! Some examples include techniques based on the evergreen Nakajima-Zwanzig 
equation from late the 1950s[17J, lla, ISSflj the hierarchica l eq uations of motion (HEOM) which provide argubly the most successful and 
robust method for calculating RpD in OQSs to date 28|, [29j , and the very recent NMQSD-ZOFE quantum master equation [26j. 



For these reasons, we will present analytic expressions for the DIF for this form of a(t) as an example in section It Table Q] represents 
a(t) in the form of equation [TJ with the pj and ttj parameters expressed analytically in terms of the inverse temperature f3 and the 
parameters of J(w), for three widely used and versatile classes of spectral distributions, each of which can essentially represent any 
physically relevant spectral distribution. The DIF for these spectral distribution functions can then be obtained in analytic form by 
substituting these pj and Qj parameters into our analytc expressions for the DIF for this form of a(i). 

Additionally, one of the most widely used forms for J(u>) is J(ui) ex w s e _ ^^ c ^[10(. TableUalso presents an analytic expression for 
a(t) for this form of spectral distribution, for certain values of s when q = 1. The DIF is then obtained in analytic form and presented 
in the supplementary material. Likewise, the DIF can be obtained in analytic form for various other spectral distribution forms in the 
same way when an analytic form for a(t) exists, whether naturally or by a numerical fit. 

The open source Mathematica program that supplements this paper takes in the size of the time step At in the discretization of 
the Feynman paths, and any form of a(t) as input; and outputs an analytic form for the DIF if possible, and quickly calculates the 
DIF numerically otherwise. This program, and the open source MATLAB program that also supplements this paper, can both readily 
calculate the bath response functions and reorganization energies of all four of the general forms of J(u) presented in table HI The 
MATLAB program can also take in any J(u>), and provide a(t) in the form of equation Q] by outputting the fitted values of {pj,ttj}. 



1 The influence functional is a functional of the Feynman paths and a regular function of all other independent variables. When the Feynman paths s(t) 
are discretized with respect to time, they can be represented by a set of variables that are constant with respect to time: {s^}. A discretized influence 
functional can therefore be thought of as a functional of discrete paths, or simply a regular function of a set of variables, and it is with respect to these 
variables that the DIF is analytic. 



RESULTS 

Setting 

The most popular OQS model is currently the Feynman- Vernon model. In the Feynman- Vernon model, the OQS s is coupled linearly 
to a set of quantum harmonic oscillators Qk '■ 

H = HoQS + -^OQS-bath + Hbath (2) 

= #oqs + J2 C « S ^ + J2 {hm*Q K 2 + ¥h^lQl) ■ (3) 

In most models the quantum harmonic oscillators (QHOs) span a continuous spectrum of frequencies uj k and the strength of the 
coupling between the QHO of frequency u) and the OQS is given by the spectral distribution function J{ui): 

j( U ) = ij2— ■*(«-««)■ (4) 

2 ^ m K uj K 
For the hamiltonian of the Feynman- Vernon model, the bath response function a(t) is the following integral transform of J(u>): 

a{i) = — I J{oS) I coth I 1 cos(wi) — i sin(wi) ) duj (5) 

1 roo J( w ) ex p(^ 



-iujt 



T J-oo sinh 



^ dw, J(-oj) = J{u) (6) 

(—J 

'dw, J(-w) = J(w), (7) 



J(w) 5 



Ti" J-oo 1 - exp(-/3uh) 
where, equation [7] can be written in terms of the Bose- Einstein distribution function with x = —j3u)h: 

jBosc-Einstcin^) = —— . (8) 

1 — exp(a;) 

In the Feynman integral formalism for expressing the RpD given this hamiltonian, all the information about the influence of 
-HoQS-bath + -^bath on the RpD of the OQS is completely captured by the (gaussian) influence phase functional (3, |j|: 

$[*(«), s'(t)] =- f f (s(f) - s'(t')) (a(f - t")s(t") - a*(t' - t")s'{t")) dt"dt' (9) 

J0J0 

= - ff $[s(t'), s '(t'),s(t"),s'(t")]at', (io) 

JoJo 

which defines the (gaussian) influence functional: 

F[s(i), S '(i)]=cxp($) . (11) 

In order to numerically calculate a Feynman integral for RpD, we discretize the Feynman paths: {s(i), s'(t)} — > {s^, s^ }£L , such that 
each s k is constant with respect to time. 

Trotter Splitting 

The simplest way to discretize the Feynman paths is to split them into intervals of equal duration, known as a Trotter splitting. This 
allows us to decompose $[s(i), s'(t)] into constituent double integrals: 



N fe -! p(k+l)At / Ak'+l)At ft' \ 

k=0k'=0 JkAt \Jk'At JkAt / 



N fc-1 r (k+l)At / r (k'+l)At _ 

(12) 

I 1„ \ J. \ I f_/ A J. I l. \ j. I 

N fc-1 

5Z 51 (4 - s k)(Vkk'st, - »?fc fc Sfc/) , (13) 

fe=0fc'=0 



where, 



Vkk' 



Vkk 



(fe+l)At r (k' + l)At 



a{t' - t")dt"dt' 
/ / a(t' - t")dt"dt' 

JkAt JkAt 



kAt Jk'At 

(k+l)At ft' 



(14) 
(15) 



At this stage, most traditionally [12, Il4, Il5|. equation [5] is inserted into the integrand of equations H4l and [T51 yielding (when 
J(-w) = -J{u))): 



_ JJ^^^I sill2( .A V2)e --A t ( fc - fc ') dw , o < W < k < N , an d 

7r J_ 00 cj 2 sinh(/3 ft "/2) 



1 
2ir 



J{uj) exp(/3^/2) 
w 2 sinh(/ 3Sw /2) 



(1-e 



-iwAt 



Or alternatively |32| . one can insert equation [5j yielding: 



iAiw)dw , < k < N . 



(16) 
(17) 



Vkk' = 2 



Vkk 



UJ 2 

J(w) 



(1 - cos(Afcj)) (cos(Aiw(fc - fc')) coth ( — — J - isin(Atcj(fc - k')))du> , < fc' < k < N , and 



(l — cos(Aia;)) coth 



i(sin(A<a;) - Attj)duj ,0<k<N. 



(18) 
(19) 



However, if a(t) = J^ p^e ° j * we can calculate equation [T4l analytically to get: 



A" 



Vkk- =4^ gsinh 2 (^At/ 2 ) e %(fe-fe')At 5 < fc' < A: < iV 



%fc 



2 Z! % ( sinh(n^*/2)e njAt/2 - -OjAtj , < fc < N 






(20) 
(21) 



Figure 1 compares the versions of Vkk' hi equations 1161 181 and [2U1 as a function of Afc = k — k' with Ai = O.OOlps , along with the 
CPU time in Mathematica, for the spectral distribution J(w) = 0.027ps 2 w 3 e^ _ "/ 2 2ps1 ' , which comes from a recent experimental study 
of quantum dots[22[ and was used in at least two recent computational studies[2| ll6(. The bath response function for this spectral 
distribution was fitted to the form of equation [T] with K ~ 4, and the fitted parameters {pj,Qj} are given in the supplementary 
material, along with a figure which demonstrates that there are no discernable differences between the least-squares fit and the original 
bath response function obtained by numerical integration. Figure 1 shows that our analytic form for rjkk' agrees essentially exactly 
with the other two forms when they are numerically integrated, but take much less time to computqj. All other analytic ij coefficients 
in this paper have been compared with essentially exact agreement against previous versions calculated by numerical integration, and 
the results of these comparisons are presented in the supplementary material. 



Strang Splitting 

A more sophisticated discretization of the Feynman paths can be made by a Strang splitting, which allows one to calculate the RpD 
much more accurately for a given number of time steps in the discretizaton of the Feynman paths. When this splitting is used, the 
above formulas for rj kk i and Vkk remain the same for < k' < k < N, but different formulas have to be used when either k or fc', or 
both, are or TV: 



2 The CPU times shown in the figure were measured on a Toshiba Satellite L300 laptop with an Intel Core 2 Duo CPU T5800 processor with a 2.00 GHz 
clock rate. 



Figure 1. Equations 1161 [181 and 1201 are equivalent, but equation [20] can be calculated much more efficiently than the other two. Since the values of 
these curves on the ordinate axis depend on the size of At (which is a parameter used for numerical computation of the Feynman integral and not 
a phsyical characteristic of the system), all of these curves are scaled by the real part of the funciton when Afc = in order to have the ordinate 
axis equal to 1 at most. Here J(w) = 0.027ps 2 o; 3 e ( ' _ "' 2 2pB ' , which is often used to study quantum dots [21, Hq, |23|. and j3 — — !-^r with T — 77K. 
Equation [TrJl was originally presented in [lj and equation [T51 was originally presented in [32|]. 
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(22) 
(23) 
(24) 
(25) 
(26) 



Again, most traditionally [12j, |l4j, [l5| equation [S] is inserted into equations |2"21 to I21H . which yields (when J(— uj) = —J(uj)): 



2 f°° JUS) exp(/3^/ 2 ) . 2 , A . N - (t AtM , 
Vno = - / ' . *\ Hr ,\ sin 2 ("At/ 4 ) e -'"(*- At /2) d ^ 



"" J-oo w 2 sinh(/ 3fiw /2) 



J{uj) exp(/ 3Rw /2) 



((1 - e - iuJ ^ 2 ) - *Ato,/ 2 ^ duj , < k < N 



^A^— sinh^) 1 
*» = « J_ ^Snf2| S ^^ M«^)e-^-^ ,0<k<N, 

VNk = 2 ~r ^^^-sin(^ V4)sin( ^ t/2)e -^( t -^-At A)dw Q k N 
ir J _ OQ to 2 sinh(^ Bw /2) 



(27) 
(28) 
(29) 
(30) 



In the case of a(t) = ^ ■ pje njt we again find analytical expressions: 



K 



j= i j 



mo = 4^ n- e W-*V*) (sinh 2 (^A t/4) ) 



?7oo = Ww 



M 



0=1 i 



%o = 4^ % 8inh(0iAt/ a ) sinh (%At/ 4 ) e %(fcAt-A«/4) 
Wfc = 4^ % sinh(^At/ 2 ) sinh (%A t / 4 ) e %(t-fcAt-Ay 4 ) 



(31) 
(32) 
(33) 
(34) 



Quasi-adiabatic displacement 



When the bath is nearly adiabatic, it is helpful to rewrite equation [3] as |ll| : 



H — HqQS — -^displacement + -ffoQS-bath + -^bath + -^displacement 



H, 



OQS 



E 



2m K ui 



Y, C - S Q + J2 (i m «^« + *m*«»lQl) +Er" 



H, 



OQS, displaced 



^c K sQ + H ] 



bath, displaced- 



els 2 



2m K w; 



(35) 
(36) 

(37) 



^displacement is called the "counter term", and can also be represented in terms of the spectral distribution by recognizing that when 
the QHOs span a continuous spectrum of frequencies lo k , we have the relation (remembering equation 2|) : 






2m K ul 7T Jo 



J(w) 



duj . 



(38) 



The influence of - ffoa s-hath + -^bath, displaced on the RpD of the displaced OQS is now completely captured by the modified influence 
phase functional 12|, [l4j, [l5|: 



$QVAI>l{s(t),s'(t)] 



' ( (s(t>) - s'(t')) («(f , t")s(t") a*(t>, t")s'(t")) + ±T [°° 
0J0 \ nn Jo Jo 



J(u) 



doj(s + (t') 2 - s-(t') 2 \)dt"dt' 



(39) 



where QUAPI stands for Quasi- Adiabatic Propagator Feynmai|f| Integral, which is the name for the resulting Feynman integral for 
the RpD of the displaced OQS when this new influence phase functional is used. When the bath is nearly adiabatic, the QUAPI is 
more accurate than a Feynman integral that does not make use of this modified influence phase functional, for a given size of the time 
step used in the discretizaton of the Feynman paths. The new term in the influence phase functional adds the following term to the 
Vkk' coefficients when k = k! : 



Aquapi 



iAi 

[At 

hit 



J(w) 



dw 



(40) 

(41) 



where in the last line we have defined the "bath reorganization energy" by A. In equations [j"6l to [T9l and l27l to l30l equation l40l can 
simply be incorporated into the integrals over ui, but our versions (equations [20l [2TJ and IBTI to l34|) are no longer analytic unless an 
analytic expression exists for A. Fortunately an analytic expression exists for A for most forms of J(w). Table U presents analytic 
expressions for some of the most popular forms of J(ui). 



3 The term 'path integral' is used more commonly than 'Feynamn integral' here, but this term is ambiguous. Currently, the first result on the search engine 
at www.google.com, when the search query 'path integral' is entered, is a Wikipedia page that currently links to three different meanings of the word 'path 
integral': (1) line integral, (2) functional integration, and (3) path integral formulation. Only the third of these is unambiguously the Feynman integral 
discussed in this paper. The 'line integral' is an integral over a path, rather than over a set of paths; and the term 'functional integral' can refer to at least 
three types of functional integrals: (1) the Wiener integral, (2) the Levy integral, and (3) the Feynman integral. 



New class of spectral distributions whose bath response functions are analytically expressed as a(i) = X/j Pi e jt - 

This formula is good for testing the spectral distribution corresponding to a bath response function obtained by a least-squares fit, 
or obtained by a truncation of the (exact) infinite series expression for a(t) table (see table U for some examples). 

J( w ) = »I(l- e -^)f;^ r (42) 

J J 

CONCLUSIONS 

We have presented general formulas for the r\ coefficients in the discretized influence functional of the Feynman- Vernon model, in 
terms of the bath response function. These formulas can then be used to derive the form of the rj coefficients presented by Makarov 
and Makri in 1995 [l4ll5| . and the form presented by Vagov et al. in 2011J32J, both which involve improper integrals of the spectral 
distribution and have at least one infinite integration bound. However, when the bath response function can be expressed analytically 
(which is almost always the case, as shown in tableUJ), our general formulas for the n coefficients can almost always be derived analytically 
too. Figure 1 shows that the evaluation of these analytic formulas can be more than 900 times faster than obtaining the exact same 
results by numerically integrating forms for the rj coefficients in the traditional ways. 

If there is a case where the bath response function can be expressed analytically, but the integrals in equation [TJJ [TS] and/or [221 to 
[25] cannot be evaluated analytically (we do not know of such a case, but cannot rule out the possibility that it exists), we recommend 
to numerically integrate these equations, rather than to numerically integrate the more traditional forms presented in equations Qj5] 
to [JjJ] and [23 to [301 This is because these integrals have finite bounds, which are very small (about the size of the time steps in the 
discretization of the Feynman paths), while the more traditional forms involve numerically integrating the spectral distribution over 
infinite bounds. 

When the bath response function cannot be expressed analytically, it may still be more convenient to obtain it by numerically 
integrating equation 0EI or [7] and then to fit it to a sum of complex- weighted complex exponentials. This is because Feynman integrals 
for calculating the reduced density operator dynamics in open quantum systems are usually used for benchmarking other techniques, 
and many popular techniques require this form for the bath response function, and when benchmarking it is always best to keep as 
many parameters of the model consistent across all methods being benchmarked. When fitting a bath response function to this form, 
equation [42] is useful for testing the quality of the fit in terms of how well the fitted bath response function corresponds to the original 
desired spectral distribution. This formula is equally useful in cases when the bath response function is naturally represented in the 
form of equation Q] (as in the examples in table H]), since it can give an indication of how many terms in the infinite series are required in 
order to be sure that one is using a bath response function which corresponds to a spectral distribution similar enough to the original 
desired one. 

SPECIAL CASES (APPENDIX) 

The results in equations [201 [23] and [3TJ to [3J] can readily be applied for many of the spectral distributions listed in table [H by 
substituting the appropriate values of pj and Qj . 

One of the most popular spectral distribution forms is currently the Lorentz-Drude/Debye (LDD) form which has been used exten- 
sively in experimental determinations of spectral distributions 24L 12a. l34| . and in various OQS studies ([SH and various other studies of 



this same system that compare to this result, including ones involving DIFs for Feynman integrals [191. I20J): 



uj A7 

■K 7 2 + LO 2 



^) = r-^-p (43) 



The first spectral distribution form presented in table [J is a generalization of this LDD form, which allows for multiple peaks to be 
included at various locations defined by {%}. In the case where only one term in the sum exists, and u) = 0, equation [53] is recovered. 
The bath response function for equation [53J can be obtained analytically in the form of equation [1] by a simple application of residue 
calculus, leading to the well-known Matsubara series. Since the Matsubara series converges extremely slowly, we have presented a series 
based on the very recent [N-l/N] Pade decomposition technique in [7j (a more clever application of residue calculus), which is also 
of the form in equation [1] but converges with much fewer terms. The label [N-l/N] denotes that the Pade approximant is a rational 
funciton where the numerator is a polynomial of degree N-l, and the denominator is a polynomial of degree N. N effectively becomes 
the number of exponential terms in a(t) that specifically arise due to this technique for expressing a(t) in the form of equation Q] (more 
details about this technique can be found in [y, [Zj)o This representation of a{i) becomes more and more accurate as N is 



N is therefore different from N that was used throughout this paper to denote the number of time steps used in the discretization of the Feynman paths. 



increased. The Pade parameters {£„,S„} are given in table [TT] The total number of terms in a(t) is given by K = 2J + N, where J 
is the number lorentzian terms of each type in J {to). 

The second spectral distribution form presented in the table below is used in studies where the J(uj) for a system is predicted 
using a molecular dynamics simulation [lj, |2l| . It is the same as the first spectral distribution form in table HI but with the prefactor 
ui replaced by tanh(^) in an attempt to remove the temperature dependence of J(u>) which is inherent due to the nature of the 
molecular dynamics technique that is used to calculate it. As for the first J(oj) presented, the bath response function for this J(uj) was 
put into the form of equation [1] by the Pade decomposition, except that in this case only the imaginary component requires a Pade 
decomposition treatment. The total number of terms in a(t) is again given by K — 2 J + N. 

The third spectral distribution form presented is known as the Meier- Tanor decomposition and has been widely used since its 
introduction in 1999 [17|. It is similar to forms used for the spectral distribution function in |27J and [3l|- Here we do not use a Pade 
decomposition treatment, but rather present the same representation of a(t) in the form of equation [1] as was originally presented by 



Meier and Tanor (using a Matsubara-style decomposition) in [17|. The total number of terms in a(t) is again given by K = 2 J + N, 



where N now denotes the number of terms contributed by the Matsubara-style decomposition. Once again, this representation of a(t) 
beocmes more and more accurate as N is increased. 

Finally, formulas [141 [15] and [22] to [26] can readily be applied to other spectral distribution functions when an analytic expression for 
a(t) exists. One very popular class of spectral distribution functions for which a(t) can be calculated analytically (resulting in a form 
different from equation [I} is the class of functions of the form J(w) oc w 8 e~( u / Uc ' . An analytic form for a(t) and the reorganization 
energy A for this class of spectral distribution functions is presented in table [I] below. The result of evaluating this formula is compared 
to the a(t) obtained by numerical integration with essentially exact agreement, and this comparison is presented in the supplementary 
material. 



Table I. Analytic formulas for the bath response functions a(t) and reorganization energies A of various selected spectral distribution functions. 
For the first three cases, K — 23 + N, where N is the number of terms in a(t) arising from the Pade (for the first two forms of J{to)) or 
Matsubara-style (for the third form of J(cj)) decomposition. All formulas presented here can be computed easily for any value of K in our open 
source MATLAB program that supplements this paper. 



A jTj 



J(u) ~"2~2j \ ^f +(u 



A jTj 



«W = Efpi^ i 



^-Tfjf ; + i_ 2 



pj= \ 



2H„n 2 
i/2-n 2 



+ i- 



J I 1 _ Y^ 
I 2_m=l 

_ V J 4A m7m 5j^j(|n m | 2 -f|) 

i->m=\ p |„ 2 -n 2 I 2 

I j ml 

r _ 7j + i^ , j e [o, j] 
- 7j - iwj , je[J + i, 2 j] 

[-l/j-, jG[2J + l,K"] 



i 6 [o, ^] 
j e [J + l, 2J] 



(•)./ +J 



2-n-gn 



A = £?A 



3 n 3 



j{„,S„) are the [N-l/N] Pade parameters for the Bose-Einstein distribution (see table Hit 



J(w) = i tanh (£*) £, J f a * i7J . , a + a J J '?. ,, ) 



*(*)=£fPie f 



J(oj) introduced by Damjanovic et aZ. [lj,[v] 



Pj 



\j . 2Aj ^jy B ti^j 
2 T 1 ,3 2^n=l ^-Q 2 

2 " r ,8 ^n=l „ 2 _fi 2 



4- j 7mA m (|n m | 2 -^ 2 ) 

i a ? ,^ — i i r 5 77t To 



It mi 



j 6 [0, J] 

j G [ J + 1, 2 J] 

ie[2J + i,JC] 



-7j + iuij , j G [0, J] 

- 7j - iu>j , j G [J + 1, 2J] 

-I/,-, je[2J + i,iC] 



l'r, 



/3 



{§ n ,2 n } =are the [N-l/N] Pade parameters for the Fermi-Dirac distribution (see table HI) 






«(*) = EfPie r 



Meier- Tanor decomposition []j 



3i-(coth(4(^+i 7j ))+i), 



^-i-( coth(|(Wj -i 7 j)) -i), 

"jA-n 



_ vr v^-^ 

S 2_, n =i 



(7 2 +(i^+<i, 1 ) 2 )(7 2 +(^j-<i, 1 ) 2 ) 



- 7j + icDj , j G [0, J] 

-7j - iu)j , j G [J + 1, 2J] 

-i/j, je[2j + i,K] 



j 6 [0, J] 

je [J + 1.2J] 

, j G [2J + 1, K] 



(■)./+, =(-h 



8 7 (7 2 



aj £f Ai 



J(w) = Au'e"''/^ 



(o 



^ C t-H) s + 1 



a(t) = -s!a;= +1 (2 + (-l)^ 1 **) 

2!R (^( 1+1^=1)) 



See |ldl | for a detailed use of this J(u>) 



) 



sHP"c) s + 1 



( (s+1) arctan(uj c t)J 
(l+o; 2 t 2 )( s + 1 )/ 2 



A = Aj=r(s) 



J(oj) = Au*e-( u / u ') q s>0,o; c >0,<j>0 



a= -w?r(^) 
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Table II. [N-l/N] Pade parameters for the Bose-Einstein and Fermi-Dirac distribution functions (first presented in [7fl and in more detail in [g). 
{H n ,£ n } can be calculated easily for all n, for arbitrary values of N in our open source MATLAB program that supplements this paper. The 
matrix A is a 2Nx2N, and A is a 2N-lx2N-l. The indices n run from 1 to half the number of non-zero eigenvalues of the corresponding matrix. 
The numbi i mvalues of these particular matrices come in pairs, for example (£„ , — £ n ); 

and for ma1 >dd numb alue will always be 0. 

{±£n} = eigenvalues^ -1 ) {±C™} — eigenvalues^ -1 ) 



A cd = 
Ka = 



Bose-Einstein 

S a,d±l 



2 x /(2c+l)(2d+l) 

s c,d±l 
2^(2c+3)(2d+3) 



S„ = (A 2 + \N) 



iC^tC* 



,ttl 






Fermi-Dirac 



A rd = 



*c,<i±l 



2^/(2c-l)(2d-l) 

S c,d±l 

2- v /(2c+l)(2d+l) 



H„ = (A 2 + iA)^ 



rc=i(C-o 
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